set.seed(5168)
## responsibility attribution ##

	library(foreign)
	library(mnormt)
	library(ordinal)

  data <- read.dta('Netherlands (2012).dta')

	model <- clmm(as.factor(Resp_attribution) ~ Prime_minister + Cabinet + Opposition + Perc_seats + Prime_minister * Perc_seats + Cabinet * Perc_seats + Opposition*Perc_seats + (1 | Party) + (1 | pid), data = data, HESS = TRUE, link = 'probit')
	summary(model)

	pars <- as.matrix(model$coefficients)
	vce <- vcov(model)
	vce <- vce[-c(12:13), ]      
	vce <- vce[, -c(12:13)]      
	
	simCoef <- rmnorm(1000, as.vector(pars),as.matrix(vce))

	summary(model)
	summary(simCoef)

	pmHi <- c()
	pmLo <- c()
	pmMed <- c()

	partHi <- c()
	partLo <- c()
	partMed <- c()

	oppHi <- c()
	oppLo <- c()
	oppMed <- c()

  otherHi <- c()
  otherLo <- c()
  otherMed <- c()

  part_pred <- c()
  opp_pred <- c()

  marg_effect <- c()
  margeff_Hi <- c()
  margeff_Lo <- c()
  margeff_Med <- c()


	for(i in 0:25){

	## partner party 
	agg <- simCoef[, 5] * 0 + 
	  simCoef[, 6] * 1 + 
	  simCoef[, 7] * 0 + 
	  simCoef[, 8] * i + 
	  simCoef[, 9] * 0 + 
	  simCoef[, 10] * i + 
	  simCoef[, 11] * 0 
		
		pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
		pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
		pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
		pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
		pred5 <- 1 - pred1 - pred2 - pred3 - pred4
		
		part_pred <- pred5 * 5 + 
			pred4 * 4 + 
			pred3 * 3 + 
			pred2 * 2 + 
			pred1

		partHi[i + 1] <- quantile(pred, 0.975)
		partLo[i + 1] <- quantile(pred, 0.025)
		partMed[i + 1] <- quantile(pred, 0.5)

	## opp party 
	agg <- simCoef[, 5] * 0 + 
	  simCoef[, 6] * 0 + 
	  simCoef[, 7] * 1 + 
	  simCoef[, 8] * i + 
	  simCoef[, 9] * 0 + 
	  simCoef[, 10] * 0 + 
	  simCoef[, 11] * i 
		
		pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
		pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
		pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
		pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
		pred5 <- 1 - pred1 - pred2 - pred3 - pred4
		
		opp_pred <- pred5 * 5 + 
			pred4 * 4 + 
			pred3 * 3 + 
			pred2 * 2 + 
			pred1

		oppHi[i + 1] <- quantile(pred, 0.975)
		oppLo[i + 1] <- quantile(pred, 0.025)
		oppMed[i + 1] <- quantile(pred, 0.5)



	#get marginal effect 
	marg_effect = part_pred - opp_pred
	margeff_Hi[i + 1] = quantile(marg_effect, 0.975)
	margeff_Lo[i + 1] = quantile(marg_effect, 0.025)
	margeff_Med[i + 1] = quantile(marg_effect, 0.5)
	
	}  
  
#seats sequence
seats <- seq(0, 25, 1)

#Now create the figure
par(mar=c(3.5,6.2,1.5,1.5))     
plot(1, 1, type = 'n',
     xlim = c(0, 25),
     ylim = c(0, 1),
     xlab = '',
     ylab = '',
     main = '', 
     axes = FALSE)
axis(1,at=seq(0,25,5),labels=T)         
axis(2,at=seq(0,1,0.50),labels=T)
polygon(c(seats, rev(seats)), c(margeff_Hi, rev(margeff_Lo)), col = rgb(red = 0.1, 0.1, 0, 0.2), border = FALSE)
lines(seats, margeff_Med, col = rgb(0.1, 0.1, 0, 1), lwd = 3)
abline(h=0,col="red")

mtext("Change in \nresponsibility attribution", side=2, line=2.5, cex=1.2)         #side = 2 means it will be on vertical axis. Line = 3 is distance from y-axis. cex is font size. 
mtext("Perceived seat %", side=1, line=2.2, cex=1.2)                   #side = 1 means it will be on vertical axis. Line = 2.2 is distance from y-axis. cex is font size. 

dev.copy(png,'Figure 2b NL 2012 - Marginal effects of party types over seat share.png', height = 395, width = 517)  
dev.off()  
